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We analyze the effect of pairing on particle transport in time-dependent theories based on the 
Hartree-Fock-Bogoliubov (HFB) or BCS approximations. The equations of motion for the HFB 
density matrices are unique and the theory respects the usual conservation laws defined by com- 
mutators of the conserved quantity with the Hamiltonian. In contrast, the theories based on the 
BCS approximation are more problematic. In the usual formulation of TDHF-I-BCS, the equation of 
continuity is violated and one sees unphysical oscillations in particle densities. This can be amelio- 
rated by freezing the occupation numbers during the evolution in TDHF-I-BCS, but there are other 
problems with the BCS that make it doubtful for reaction dynamics. We also compare different 
numerical implementations of the time-dependent HFB equations. The equations of motion for the 
U and V Bogoliubov transformations are not unique, but it appears that the usual formulation is 
also the most efficient. Finally, we compare the time-dependent HFB solutions with numerically 
exact solutions of the two-particle Schrodinger equation. Depending on the treatment of the ini- 
tial state, the HFB dynamics produces a particle emission rate at short times similar to that of 
the Schrodinger equation. At long times, the total particle emission can be quite different, due to 
inherent mean-field approximation of the HFB theory. 
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I. INTRODUCTION 



Pairing is essential to the global description of nu- 
clear ground-state and low excited state properties; the 
Hartrec-Fock-Bogoliubov (HFB) and Hartrce-Fock aug- 
mented by BCS (HF-I-BCS) theories are in common use 
to treat the pairing degrees of freedom Also m nu- 
clear reactions, many phenomena are expected to be in- 
fluenced by pairing correlations: collective motion, fu- 
sion, fission, transfer reactions and nuclear break-up. 
The obvious candidate theory to treat these effects is 
the Time-Dependent Hartree-Fock Bogoliubov (TDHFB) 
theory 0, and there has been much effort in the last 
decade to apply it. However, the TDHFB theory turns 
out to be much more complicated to implement than 
the corresponding time-dependent Hartrec-Fock theory, 
and the applications have been mainly performed in its 
small amplitude limit, the Quasi-particle RPA (QRPA) 
[sl-fioj. However, most of the phenomena quoted above 
are far from small amplitude excursions from the ground 
state, and transport theories able to treat Large Am- 
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plitude Collective Motion (LACM) arc mandatory. Re- 
cently, several groups have applied the TDHFB pl| - [l3j 
to nuclear dynamics. An approximate version of theory, 
called TDHF-fBCS, has also been considered [H. We 
will discuss its properties as well. Most recent applica- 
tions have been to small-amplitude collective motion in 
nuclei, where the theory is equivalent to the quasiparticle 
random-phase approximation (QRPA). Still, it is impor- 
tant to understand and solve the theory here as a first 
step towards treating large-amplitude motion. 
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The aim of the present article is first to present from 
a rather general point of view different dynamical theo- 
ries that incorporate pairing correlations. We find that 
the TDHFB has many good properties, but it is difficult 
to find further simplifying approximations. We find that 
the TDHF-I-BCS approximation leads to a break-down 
of the continuity equation. This failure might induce se- 
rious difficulties in the description of physical processes. 
The second aim of the paper is to test various implemen- 
tations of the pair theory in a model problem. For this 
purpose, we examine a one-dimensional model of parti- 
cle evaporation for which we also have a numerical exact 
solution. 
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II. FORMALISM 

The TDHFB theory and the TDHF+BCS approxima- 
tions have been recently appHed to nuclear physics in 
Ref. [TT| - [l^ . In this section we will briefly summarize 
the main features. 



A. The TDHFB theory 

The main equations of TDHFB theory, Eq. be- 
low, can be derived in at least two different ways. One 
way is from the general variational principle. 



{^{t)\ihdt - H\'9{t))dt, 



(1) 



see e.g. Ref. [l5[. Here H denotes the Hamiltonian 
and l^*) is a HFB wave function. The equations may 
also be derived by demanding that the operators for the 
ordinary and anomalous densities, pij = a^a^ and kij = 
ajUi respectively, satisfy Ehrenfest's theorem: 

ihdt{^\a]a,\^) = m[a]a,,H]\^), (2) 
ihdt{^\a,a,\^) = {^\[a,a„H]\^). (3) 

In the following, we will further assume that _ff is a two- 
body hamiltonian. 



(4) 



ijkl 



where v denotes the anti-symmetric two-body matrix el- 
ements. The derived TDHFB equations are: 



ih — p = hp — ph + kA* — Ak* , 
at 



(5) 



ih—K = hn + Kh* + A{1 - p*) ~ pA. (6) 
at 

Here p, k, h and A arc all matrices of dimension equal to 
that of the single-particle space. The matrices h and A 
are the mean-field and pairing field of the Hamiltonian, 
defined as 



kl 



= \ ^VijklKkl- (7) 



kl 



The dynamical equation can be recast in a more com- 
pact form by introducing the generalized density matrix 
TZ and generalized single-particle hamiltonian %: 



n = 



p K 
-K* 1- p* ' ' 



n = 



h A 

-A* -h* 



(8) 



With these definitions the equation of motion be- 
comes [m, Eq. 9.61a] 



7^= [H,TZ] 



(9) 



This generalizes the usual TDHF picture by replacing 
the one-body density by 7Z. Similarly to the TDHF case, 
the generalized density satisfies Ti? = TZ and has only 
eigenvalues equal to zero and one [l^, [l^ . 

Going back to ordinary Hartree-Fock theory, it is com- 
putational advantageous to factorize the density matrix 
and express it as a sum over the contributions from oc- 
cupied orbitals to obtain equations of motion for the in- 
dividual orbitals. There is no obvious advantage for the 
factorization in TDHFB because all of the single-particle 
orbitals in Fock space contribute to the generalized den- 
sity matrix TZ. Nevertheless, the factorization is usually 
applied to obtain the actual equations to be solved nu- 
merically. To write equations in this form, one needs an 
explicit form of the Bogoliubov transformation. 



(10) 



The density matrices are expressed as p = V*V'^ and 
K — V*U'^ , and the generalized density matrix is 



n = 



V* 

u* 



(11) 



One can then easily see that Eq. ([9]) will be satisfied if 
we require the {[/, V} matrix be a solution of 



iTi 



h A 

-A* -h* 



(12) 



The numerical solution of the TDHFB equations are usu- 
ally carried out in this representation [l3. fl3j. However, 
it should be remembered that there are redundant vari- 
ables in the {[/, V} representation corresponding to uni- 
tary transformations of the quasiparticle basis, and in 
fact Eqs. are not unique. 

We may derive another form of the TDHFB equa- 
tions as follows. The wave function j^*) at any time t 
is the quasi-particle vacuum associated with the Bogoli- 
ubov transformation that transforms the physical vac- 
uum to |5'(t)). In that representation, the Hamiltonian 
has zero-, two-, and four-quasiparticlc terms that can act 
on \'^{t)) [l^. Neglecting the four-quasiparticlc excita- 
tion amplitudes, the result is 



H\^{t)) H'{t)\^{t)) 



{H) + WH^$PHt)pl{t) 



afi 



l*W>(13) 



with [il, Eq. (E.22)] 

According to the Thouless theorem, any state of the 
form (1 + Za/3/3l^(3\)\^) can be expressed as a new 
quasiparticle vacuum [l6|. To lowest order in Z, the 
Bogoliubov transformation to the new vacuum from the 
physical vacuum is given by 



(U' V*' ) ^ {U V* ) 



1 Z* 
Z* 1 



(15) 
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Wc can thus derive an equation of motion by demanding 
that the changes in U, V just match the two quasiparticle 
excitations generated by H. After a lengthy but straight- 
forward derivation, it may be shown that the correspond- 
ing equations of motion for U and V can be written as: 

r ihdtU = pK^U -Kh*V 

- kA*U + pAV 

(16) 

ihdtV = -{1- p*)h*V - K*h'<U 

- K*AV- {1- p*)A*U 

These equations differ from (|12|) but nevertheless lead to 
the same TDHFB equation for the generalized density. 



time-dependent orbitals. For computational reasons the 
phase is removed by integrating 

iTidtln) = {h[p] - Vk)\n) (21) 

instead of Eq. ([TH]), with ijkit) = {ipk{t)\hHF\^k{t)) . At 
the same time, Eq. (pOj) is replaced by 

i?i^Kfe = Kfc(?7fc + ?7j:) + Afe(l - 2nfc). (22) 

Finally, we mention that the TDHF-I-BCS approximation 
was found to work well with a separa ble p airing interac- 
tion and in the small amplitude limit |14l |. 

C. Conservation laws and equation of continuity 



B. The TDHF+BCS approximation 

The TDHF+BCS treatment of pairing dynamics is mo- 
tivated by the simple form the wave function has in the 
BCS approximation, 



I*) 



n 

fc>0 



Uk 



Vkalal 



(17) 



The TDHF+BCS approximation may be derived from a 
variational principle [l7| or by an approximate reduction 
of the TDHFB equations [1J|. For the reduction of the 
TDHFB equations, we first note that wave function can 
be put into BCS form at any fixed time by transforming 
the U, V matrices to the canonical basis. In that basis, p 
is diagonal and k matrix is zero except for one element 
on each row (or column) representing the pair ii. Assum- 
ing that the A matrix has the same structure as k, Ref. 
[14| shows that the TDHFB time evolution preserves the 
same canonical structure with orbitals that evolve by the 
mean field Hamiltonian, 



ihdt\ipk) = h\ipk). 



(18) 



where h has been defined in Eq. ([7]). The equations of 
motion for p and k in this time-dependent basis are^ 



'i-^-^i^k = +Afc(l-2nfc). 



(19) 
(20) 



Here , and Afc are short-hand notations for pkk , f^kk 
and A^^ respectively. 

One technical point should be mentioned. When Eq. 
([T8| is integrated, there is an irrelevant phase factor 
exp(— I J* {(fkit')\hHF{t')\(pk{t'))dt') introduced into the 



^ Note that these equations slightly differ from those from [l4j . 
due to the definition of here. 



Since the TDHFB density matrix satisfies Ehrenfest's 
theorem, it is trivial to show that the conservation laws 
for one-body observables are respected by the TDHFB 
dynamics. It was also shown that conservation laws for 
important observables such as particle number are satis- 
fled in TDHF-BCS [lj|. However, for transport we are 



interested in local conservation laws as well. In particu- 
lar, if the interaction is local the coordinate-space density 
n{x,t) should satisfy the equation of continuity, 

dn{x, t) 



dt 



(23) 



where j{x) is the particle current. Assuming a local in- 
teraction, Eq. ([25)1 may be derived from Ehrenfest's the- 
orem, evaluating the commutator on the right hand side 

as 

t;2 



S/-j{x)^{[h{x),H]) 



2m 



{[h{x),V^]). (24) 



This is sufflcient to guarantee that TDHFB obeys the 
equation of continuity under the stated condition. Unfor- 
tunately, this is not true for the TDHF+BCS dynamics. 
Within TDHF+BCS, the local density is given by 



n{x, t) 



y n,;(t)|(^,(x,i)|', 



(25) 



and its evolution satisfies 
dn{x, t) 



dt 



= ^'n^{(p*{x,t)^t(f^{x,t)) + ifi{x,t)dtf*{x,t)) 



+ 



\ip^{x,t)\^dtn,{t) 



(26) 



The first terms on the left are just the evolution of the 
orbitals under a mean-field potential, and so the same 
reduction applies as in Eq. (|23p . The result is 



dn{x, t) 
dt 



= -V-J(.T,t)+^|(/,,(x,i)|2 



dni{t) 
dt 



(27) 



Thus continuity cannot be guaranteed unless the oc- 
cupation numbers are fixed. We will see below that 
TDHF+BCS can produce unphysical density oscillations 
when the occupations are allowed to vary. 
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III. APPLICATION TO PARTICLE 
EVAPORATION 



Recently two of us (DL and KW) began investigating 
the effect of pairing on particle evaporation, and obtained 
the results shown in Fig. [T] Skipping over the details, the 
number of particles escaping an initially excited nucleus 
is shown as a function of time using either the 3D-TDHF 
code of Ref. [l8l - [20| or an upgraded version including 
pairing using the TDHF+BCS theory proposed in Refs. 
M UM- As we can see, the standard mean-field cal- 
culation presents the expected long time decay due to 
particle evaporation 2l| . When pairing is included, the 
number of particles in the nucleus first decays and then 
starts to oscillate. Clearly, this result is unphysical. It 
was this unphysical result that motivated us to under- 
take the present more general study. For the present 
article, we consider a more simplified Hamiltonian that 
permits us to compare a number of approximations with 
each other and with a numerically exact solution. In the 
present article, we investigate whether the observed prob- 
lem is systematic in theories where pairing is included or 
if it comes from the specific treatment of pairing in the 
TDHF+BCS approximation using zero range interaction. 
Our study is also the occasion to benchmark different 
theories, TDHF+BCS and TDHFB, to describe particle 
emission. 
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FIG. 1: Top: Schematic illustration of neutron evaporation 
from a nucleus of O^^ excited by a monopole boost at t=0. 
Bottom: Number of neutron inside a sphere of size 10 fm 
around the nucleus as a function of time obtained with TDHF 
and TDHF+BCS (from 



A. A one-dimensional model 

For comparing the different treatments of pairing dy- 
namics, we consider a one-dimensional system composed 
of iV particles in a box with x in the range — ^max < x < 
^va&y. and a Hamiltonian of the form 

N 



H 



E 



2m 



+ U{x,) 



N(N-l)/2 
i<j 



(28) 



Here, Pcnaj denotes the spin-exchange operator. The po- 
tential U (x) is taken to be a Woods-Saxon well centered 
at the origin: 



Uix) 



l+exp[(|a;|-Xo)/a] 



(29) 



The two-body interaction v{x — x') is taken to be a finite- 
range Gaussian 



v{x — x') = Wo exp 



{x - x'f 



(30) 



In the limit where the range (Tq goes zero, v{x — x') is a 
contact interaction and our model is similar to the model 
considered in Ref. [2^ to analyse the onset of vortices 
in rotating Fermi gas using TDHFB. The advantage of a 
finite range is that it does not have to be renormalized 
for use in BCS or HFB. 

The TDHFB is formulated in a Fock space and the 
space has finite dimension in numerical implementations. 
Our particle creation and annihilation operators 
are defined on a uniform mesh of points {x} with spacing 
Ax; a =t 01' i is the spin label. Then we can write the 
quasiparticle transformation as [23j 



Pi 



\x'^{ua{x,t)i!\{x) + Va{x,t)lpi{x)^ . (31) 

X 

ia;^ {ua'{x,t)'4)\{x) + Va'{x,t)%l)^f{x)^ .(32) 



In the following, we will use the convention Ax X^a; ~^ 
X^ and not distinguish between the quasiparticle sets 
a and a', with the property Ua'{x,t) ~ Ua{x,t) and 
Va>{x,t) = ^Va{x,t). The discretized time-dependent 
equations in version Eq. (|12p of TDHFB take the ex- 
plicit form 



in — Ua[X,t) 

ot 



and 

i^-^Va{x,t) 
ot 



(2) 



2m(Aa-)2 



+ ;7(a:) + r(a;) } Ua{x,t) 



^ A(a;,a:')wa(a;',<) 



(33) 



2m{Axy 



+ Uix)+T*{x)}v^{x,t) 
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-^A*(x,a:'K(x',t). (34) 

x' 

with 

r(x) = ^v{x-x')p{x',x'), (35) 

x' 

A(x,x') = v{x — x')k{x,x'). (36) 

Here aI^'' is the second-difference operator, A'^^)0(z) = 
0(z + f)-2(/)(i) + (/)(i-f). 

The normal and anomalous density matrix are given 

by 

p{x,x) = Y,K{xM (37) 
k{x,x') = ^<(x)M„(a:')- (38) 



B. Exact solution for the two particle case 

One interesting aspect of the model considered here is 
that for two particles it can be solved exactly numerically. 
Indeed, assuming that the system is a spin singlet, the 
two-body wave-function reads: 

$(a;i,cri,a;2,(T2) = -^{^cji^^a^i - 5aii5a2^) <i}{xi,X2i?>^) 

where 4>{xi, X2) is a symmetric function that satisfies the 
Schrodinger equation: 

ih^4>{xi,X2) = [/i? + h^2+ - 2^2)] 0(a;i,X2), (40) 

Since the discussion here might be applied not only 
to nuclear systems but also to other field of physics 
like condensed matter or atomic physics, we consider 
here reduced units. The length, time-scale and energy 
scale given below are respectively written in units of Ax, 
mA.x'^ /fi and It? / [mA.x'^) where Aa; is the discretiza- 
tion mesh step. Accordingly, all quantities below will 
be presented without specific units. The parameters of 
the central potential are set to a = 2, Xq = 4.5 and 
(To = 2.5 and the initial harmonic constraint is taken as 
A = 6.173 X 10~^. Three interaction strength vq equal 
to -1.096 X 10-2, -3.344 x lO'^, and -6.280 x IQ-^ are 
considered. The three cases will be referred respectively 
to case (a), (b) and (c) below. In each case, the depth 
of the Woods-Saxon potential has been adjusted to get 
the same binding energy E = —2.2 x 10-^^ leading to 
Uo = -2.7 x 10-2, -1.929 x lO'^ and -7, 716 x lO'^ re- 
spectively. For cases (a) and (b) the interaction is below 
the strength needed for a condensate in the HFB or BCS 
theory at a mean particle number of two. 

An illustration of the two-body density matrix ob- 
tained in different level of approximation for the case 
(a) arc shown in Fig. [21 Due to the attractive nature 



of the two-body interaction used, the two-body density 
presents a clear correlation along the axis {xi + X2)/2 
that is completely neglected at the Hartree-Fock level. 
Such a correlation is partially recovered when pairing is 
included in the HFB or BCS theory. 



Exact Hartree-Fock 




Xi Xi 

FIG. 2: S — component of the two body density matrix 
p^^\xi t, a;2 4.) in 10"^ (unit of length)"^ at time t = for 
the four theories studied here. This figure correspond to the 
set of parameters (c) (see text). 



C. Some numerical aspects for dynamics with 
pairing 

It is important to integrate the time-dependent equa- 
tions of motion with a high-order method, because wave 
function conditions such as normalization and conserved 
quantities such as energy can be easy lost. The time scale 
for single-particle motion and direct reactions is several 
thousand of units of time, and we require numerical accu- 
racy up to those times. For most of the results we present 
below, we have used the fourth-order Runge-Kutta algo- 
rithm (RK4). The calculations in Ref. [Ijl on the other 
hand use a sixth-order Adam-Bashford algorithm, and 
we have tested that as well. 

Typically, we take a box of dimension X^ax ~ 500, 
giving the HFB matrices a dimension of AXmax/ = 
2000. The single-particle Hamiltonian has a range up to 
~ 2 , which requires a fairly small time step. We take 
At = 0.263. 



1. Ehrenfest vs Thouless equation of motion 

As has been stressed in section III Al the equation of 
motion on the (uq,Uq) components are not unique. We 
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have implemented two of the formulations below, namely 
the "Ehrenfest" (Eq. ^) and the "Thouless" ( Eq. 
(|16p ) equations of motion. The numerical integration 
can be carried out very accurately using each version 
of the equations. We found that the Thouless equation 
has a better precision than the Ehrenfest equation using 
RK4 at a fixed time step. However, it turns out that the 
Ehrenfest formulation is three or four times faster than 
the Thouless one, due to the smaller number of matrix 
operation in Eqs. compared to Eqs. ([TC]). Since the 
computational time is a crucial aspect of the numerical 
treatment, the standard Ehrenfest equation is a better 
choice. The density formulation (Eq. ([5]|6])) would have 
a similar number of matrix operations to the Thouless 
formulation, but we have not investigated the numerical 
performance of this third alternative. 



2. Imaginary absorbing potential 

Particle loss is monitored by computing the number of 
particles having |a:;| < Xo/2. Particles can be reflected 
from the edges of the box and obscure this measure of 
evaporation, so we have to add an absorbing potential hi 
near the edges. As mentioned in Ref. [l2|, the specific 
form of the absorbing potential is not obvious, because it 
should decrease the particle number without affecting the 
normalization of the wave function. It can be shown that 
the following prescription satisfies these requirements, 

h — phi —A — nhi 1 ( " 1 

dt[v ) ~ 1^ -A* + K*h, -h* + (1 - p*)h, )\v )■ 

In particular, the above equation preserves the unitarity 
property uu^ + v*v^ = 1. 

In applications below, the imaginary potential is taken 
as 

hi{x) = for \x\ ( {Xmax - Xim), 

, i \ -fr 1^1 ~ X,nax ^ Xim r | | \ / \a \ 
h^{x) = iVim for \x\ ) [Xmax-Xim) 

with Xmax = imax/2, Vi,n = -7.716 X 10"^ and Xir,^ = 

37.5. 

We have compared TDHFB evolution in small box in- 
cluding the imaginary potential with the corresponding 
evolution in very large box to check that the present 
method is a practical way to suppress the reflected par- 
ticles. We also found that a simplier prescription is also 
adequate for our purposes. Namely, one can apply the 
imaginary potential to the v amplitudes alone, with the 
equation of motion 

'^di{l) " (-A* -h*^^) (")• 

This prescription violates unitarity, but the results using 
it could not be distinguished from the correct evolution. 



D. Particle evaporation 

To simulate an evaporating system, we start with a 
wave function that is constrained to be largely inside the 
potential well U{x). This is achieved by adding a small 
harmonic constraining field to the Hamiltonian and 
solving for the HFB ground state. At time < > 0, the 
harmonic constraint is removed inducing a monopole os- 
cillation of the system that is eventually damped out by 
particle evaporation. This is illustrated in Fig. [31 show- 
ing snapshots of the density at different times with the 
system evolved with the TDHF equations of motion, ie. 
without pairing. 

t=0 t=656 t=1313 



0.4 



0.2 - 



0.0 ^ 




-300 300 -300 300 -300 300 

XXX 

FIG. 3: Evolution of the local one-body density n(x) of a 
system oi N = IQ particles. The system is initially confined 
in a harmonic trap. At t > 0, the external constrained is 
relaxed. 



1. Comparison between the exact solution and TDHFB 

In this section we will compare the particle emission of 
TDHFB with that given by the two-particle Schrodinger 
equation, solved numerically. One should not expect 
close agreement under all conditions for two reasons. The 
total emission probability in the final state can be calcu- 
lated easily in the Schrodinger dynamics by taking the 
overlap of the initial state with the bound solutions. The 
TDHFB dynamics on the other hand may have no bind- 
ing when the average particle number on the nucleus be- 
comes small. 

It is also not possible to set the initial conditions for 
the HFB wave function to correspond exactly to the two- 
particle wave function of the Schrodinger equation; one 
sees this already in Fig. 2. As described above, the ini- 
tial state for the Schrodinger equation is squeezed ground 
state, namely the lowest state of the two-particle system 
in the presence of a harmonic external potential. A cor- 
responding HFB wave function could be constructed by 
using the BCS form of the wave function and requiring 
that it have the same one-particle density matrix. This 
turns out to not work well, due to high momentum com- 
ponents in the wave function that are not properly con- 
trolled by the HFB pairing field. We found that a better 
prescription is to make a corresponding squeezed ground 
state in the HFB treatment. We use this prescription for 
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the comparison shown below. 

We measure the number of particles inside the system 
by the quantity 



Nit) 



2n{x, t)dx, 



(41) 



where Xbox here is taken as 100. Note that the system 
is centered at x = 0. The evolution of N{t) for several 
cases is shown in Figure |4l comparing the exact results 
with the HFB approximation. 




2000 
time 



3000 



4000 



FIG. 4: Number of particles evaporated from an initially com- 
pressed system with initially A'' = 2 as a function of time ob- 
tained with the exact (solid black line) and TDHF (filled green 
triangle) and TDHFB (open blue squares). Results with dif- 
ferent two-body interaction strengths (case (a), (b) and (c)) 
are respectively shown from top to bottom (see text). 

In the case (a) and (b), TDHF and TDHFB are iden- 
tical. Indeed, the minimization of HFB equation to get 
the initial state leads to a pure Slater determinant state. 
In these case, the TDHF evolution is very close to the ex- 
act solution. Note that, in this regime, the evaporation 
is dominated by the mean-field contribution and pairing 
has a weak effect on particle emission. As the interaction 
strength increases, the TDHFB and TDHF results starts 
to deviate from each other as well as from the exact evo- 
lution. As the interaction strength increases, the role of 
pairing and, more generally, correlations on evaporation 
becomes more important. The TDHF evolution largely 
underestimate the emission in case (c). This stems from 
the fact that mean-field is not able to properly describe 
the diffusion of the occupation probability around the 
Fermi energy in the initial state and the dynamical scat- 
tering of single-particles during the evolution induced by 
correlations beyond the Hartree-Fock. In bottom panel of 



figure m the lack of evaporation in TDHF is due to the 
fact that all initial occupied states can be decomposed 
onto bound states of the corresponding mean-field. A 
similar situation occurs for the = 10 case presented 
below. 

A precise study of the strongest coupling case (case 
(c)), which is the only case above the HFB threshold 
for the initial state, shows that the time scale associated 
to particle evaporation is properly accounted for in TD- 
HFB. This could indeed be seen in bottom part of figure 
m where we see that the time at which N{t) starts to 
decrease is the same in the exact and in the TDHFB 
case. This shows that the time-scale associated with the 
evaporation process is the same in the exact and TD- 
HFB case. In the long time limit, TDHFB overestimates 
the average number of emitted particles. Accordingly, 
it could be anticipated that the internal motion of the 
system is more damped in the latter case than in reality. 
We indeed have checked that the damping width of the 
monopole resonance is larger in TDHFB compared to the 
exact solution. 

It should be noted that the approximation leading to 
TDHFB can only be justified for the short-time evolu- 
tion. Indeed, even starting from a quasi-particle state, 
correlation beyond TDHFB might built up in time, like 
for instance four quasi-particle excitations. 



2. Comparison between the exact solution and TDHF+BCS 

Here, the results obtained by using the TDHF+BCS 
equation of motion discussed in section III Bl are pre- 
sented. In figure [5l an illustration of the result obtained 
in the case (c) is shown. 




4000 



FIG. 5: Number of particles evaporated from an initially com- 
pressed system with initially A'' = 2 as a function of time ob- 
tained with the exact (solid black line) , TDHF (dashed green 
line) and TDHF-fBCS theory (thin red line) in the case of 
parameter set (c). 

Independently of the set of parameters used in the 
model case, it is generally observed that the TDHF-I-BCS 
theory leads to unrealistically fast early emission of parti- 
cles compared to the exact case. This fast emission seems 
to be a generic feature of the BCS approach as illustrated 
in Figs, ini where A^ = 10 particles are considered. 
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Independently of the set of parameters used in the 
model case, it is generally observed that the TDHF+BCS 
theory leads to unrealistically fast early emission of parti- 
cles compared to the exact case. This fast emission seems 
to be a generic feature of the BCS approach as illustrated 
in Figs. [6] where N = 10 particles are considered. This 
observed fast emission is very likely connected the prob- 
lem of applying BCS when continuum states are present 
in the wave functions . The BCS ground state has a un- 
physical gas of particles in the continuum rather than an 
exponential decay into the vaccum This was one 

of the historical reason why HFB was preferred to BCS 
in nuclear structure studies. It is of course possible to 
reduce the continuum problem by truncating the num- 
ber of single-particle states that contribute to pairing. 
However, we do not know any systematic way to carry 
this out without reference to more reliable calculational 
methods. 

In studies dedicated to nuclear structure, this is gen- 
erally circumvented by reducing significantly the num- 
ber of single-particle states that contribute to pairing. 
Then, only states with single-particle energy within a 
given range AE around the Fermi energy are used, where 
AE is of the order of few MeV. In FiguresHEl this restric- 
tion has not been made and a large set of single-particles 
is retained. If the energy window AE is reduced, the 
time-scale associated to particle evaporation is increased 
and eventually becomes more consistent with the exact 
dynamics. Conjointly, the asymptotic number of evap- 
orated particles is significantly reduced and approaches 
the TDHF case as AE goes to zero. 
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FIG. 6: Number of particles evaporated from an initially com- 
pressed system with initially A'' = 10 as a function of time 
obtained with the TDHF (solid line), TDHFB (dotted line) 
and TDHF-I-BCS theory (dashed line). 



It should be mentioned that in realistic three- 
dimensional calculations, there is no flexibility in the se- 
lection of single-particle states contributing to the dy- 
namics. Indeed, static calculation are already made with 
a specific choice of single-particle space in such a way that 
with an effective force in the pairing channel, the gap has 
a reasonable value. Accordingly, the dynamics should be 
made with the same set of single-particles states has is 
already done in Ref. [l^ . 



3. Spurious oscillation in TDHF+BCS theory 

In the long time evolution, oscillation of the number 
of particles, similar to those displayed in figure [T] are ob- 
served in TDHF-KBCS, see Figs. [l]and[6l Such oscilla- 
tions are absent in the TDHFB theory. From the appli- 
cation presented here, we can conclude that the spurious 
oscillations is a generic effect in TDHF-I-BCS. It occurs 
even if a finite range interaction is used. Finally, this 
problem is solved when TDHFB is used. 

To better characterize the oscillation, N{t) can be ex- 
pressed in the canonical basis as 



(42) 



where rii (t) and Pi (t) denote respectively the occupation 
numbers and the probability of the canonical orbital i 
inside the box: 



\(pi{x,t)\'^dx. 



(43) 



An illustration of N{t) for the two particle case (c) is 
given in figure [T] The observed evolution is mainly due to 
the evolution of the two closest levels below and above the 
particle emission threshold labelled respectively by "1" 
and "2". These two levels verify ni(t)-|-n2(t) — 1. During 
time evolution, the unbound level is continuously emitted 
while the bound level remains in the box, i.e. Pi{t) = 1. 
Assuming, that only these two levels contribute to the 
particle emission, an estimate N'{t) of the number of 
evaporated particles is given by ^: 



N'{t) = 2[n,{t)P,{t)+n2{t)P2{t)] 



(44) 



The evolution of ni{t) and Pi{t) for i = 1,2, as well as 
N'{t) are shown in Fig. [7] attesting for the validity of the 
two-level approximation. As seen in bottom part of this 
figure, N'{t) is very close from its exact value N{t) and 
oscillations arc due to oscillations in occupation numbers 
Such oscillations of occupation numbers are expected 
in any theory beyond TDHF, including the TDHFB 
and/or exact evolution (see Fig. [8|). However, these theo- 
ries do not lead to unphysical evolution of particle num- 
ber. The difference between TDHF-t-BCS and the two 
other theories stems from the approximation made to get 
the equation of motions. Indeed, by neglecting the off- 
diagonal matrix elements of the pairing field, the single- 
particle evolution reduces to self-consistent mean-field 
dynamics, similar to the TDHF one. The effect of corre- 
lation only enters into the occupation numbers evolution, 
and only affects the single-particle evolution through the 
density dependence of the self-consistent mean-field. 



Note that the factor 2 here comes from the initial degeneracy 
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FIG. 7: Top: Evolution of occupation numbers of the two 
closest states above (state 2, dashed line) and below (state 1, 
solid line) the Fermi energy as a function of time obtained in 
TDHF+BCS (parameters set (c)). Middle: Evolution of the 
corresponding portion of the wave-function remaining inside 
the box. Bottom: Evolution of N{t) (thin line) and of N'{t) 
(thick line) as a function of time. 



Usually, correlation is expected to induce a mixing of 
single-particle states. Indeed, the evolution of the one- 
body density matrix in the presence of correlation is given 
by: 

ifi^ = [hip),p]+Tr2[v,2,Ci2], (45) 

where h(p) is the mean-field of the correlated state while 
Vi2 and Ci2 denotes the two-body interaction and corre- 
lation matrix respectively (see Ref. [25| for more details) . 
In both TDHFB and exact solution, the second term in- 
duces an extra mixing of single-particle states that is ne- 
glected in TDHF-fBCS. It turns out that this mixing is 
essential to compensate the possible oscillations in oc- 
cupation numbers. This is clearly illustrated in Fig. [5] 
where the quantity Pi (t) are shown to oscillate coherently 
with Ui (t) in the exact case (similar behavior is observed 
in TDHFB evolution). 



4- Link with the break-down of continuity equation in 
TDHF+BCS 

Starting from the expression (j27p derived for 
TDHF+BCS, the evolution of particle number inside the 



0.06 




0.0 I ■ 1 ■ 1 ■ 1 ■ -J 

500 1000 1500 2000 



FIG. 8: Top: Evolution of occupation numbers of the three 
main single-particle canonical states contributing to the par- 
ticle evaporation for the exact dynamics. The corresponding 
values of Pi{t) are shown in the bottom part. 

box is given by 

^ = - J dW{j{x,t))dx 

i ^ ^ 

Introducing two sets of real functions Ri{x, t) and Si{x, t) 
for each wave-packet such that: 

ip^{x, t) = R4x, t) cxp(iS',(x, t)/h), (47) 

and making use of partial integration technique, the first 
term in eq. (|46| can be recast as: 

j div{j{x,t))dx ^ 2'Y^ni\ipi{Xhox,t)fvi{Xhox,t), 

where Vi denotes the local velocity of the particle defined 
through Vi{x,t) = VSi{x,t)/m. 

This term is the expected physical term expected to 
appear in any well defined transport theory that relates 
the number of particles inside the box to the flow of par- 
ticles outgoing at the boundary of the box. However, due 
to the presence of the second term in eq. (|46|). oscilla- 
tion of occupation numbers that are not compensated by 
oscillation of the probability Pi{t) (see figure [7]) lead to 
spurious behavior of the particle number. The only way 
out to avoid this problem in a TDHF-I-BCS approach is 
to freeze the occupation numbers during the evolution. 

5. TDHF+BCS with frozen occupation numbers 

To incorporate pairing in a transport model we are 
facing the difficulty that the TDHFB theory is very de- 
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manding numerically. A possible solution to this diffi- 
culty, would be to use the simpler TDHF+BCS approach. 
However, in view of preceding sections, the approxima- 
tion made to obtain TDHF-I-BCS leads to unphysical be- 
havior especially when continuum plays a significant role: 
strange behavior of particle emission, gas problem. 

We have seen in section IIII D 31 that the pathologies 
of TDHF-I-BCS comes from the evolution of occupation 
that should normally be accompanied by a consistent 
mixing of the single-particle states along the dynamical 
path. This approximation does not seem to be critical in 
the study of static properties of nuclei and most often, 
for not too exotic nuclei, BCS theory provides a fairly 
good approximation to HFB. 

A simple prescription to avoid non-physical evolution 
in TDHF-fBCS is to assume that the occupation num- 
bers are frozen during the time-evolution, this approxi- 
mation is called hereafter frozen occupation approxima- 
tion (FOA) . An illustration of the FOA effect on particle 
evaporation is shown in figure [5] (dashed line). As antici- 
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FIG. 9: Evolution of the number of particles evaporated 
from an initially compressed system of A*' = 2 particles. The 
exact result (thick solid line) is compared to the TDHF-I-BCS 
with (dashed line) and without (thin solid line) the frozen 
occupation number approximation. The simulation has been 
made with the same parameters set as the lower panel of figure 

Bl 

pated, spurious oscillations of the particle number evap- 
oration disappear in the FOA. It turns out, that, for the 
specific set of parameters used in the example of figure HI 
the asymptotic number of particle evaporated is in very 
good agreement with the exact case, much better than 
the TDHFB solution (see figure S]). However, the agree- 
ment depends on the parameters that arc used and not 
systematic conclusion can be drawn. Note that this ap- 



proximation has already been used in realistic calculation 
for example to study dipole giant resonances [26j . 



IV. SUMMARY 

In this article, different transport theories able to in- 
corporate pairing arc discussed. One important conclu- 
sion is that theories like TDHF+BCS where the conti- 
nuity equation is not respected can lead to unphysical 
results. More specifically, the effect of pairing on parti- 
cle emission has been analyzed here using a simple one- 
dimensional Hamiltonian that can be solved exactly for 
the case of two particles. From the systematic study we 
have made by changing the interaction strength and/or 
particle number, pairing does affect significantly the par- 
ticle number emission. While the TDHF approach gen- 
erally underestimate significantly the number of emit- 
ted particles, an enhancement of particle evaporation is 
observed when pairing is included. This effect is auto- 
matically included in both TDHFB or TDHF+BCS the- 
ory. Only TDHFB provides a good description of particle 
emission at short time but might deviates from the exact 
dynamics at longer time due to accumulated correlation 
effects beyond this approach. 

While the asymptotic number of emitted particles is 
quite reasonable. TDHF+BCS leads to unphysical rapid 
emission and spurious oscillations of the number of emit- 
ted particles. The direct use of TDHF+BCS, that would 
be highly desirable from the practical point of view, 
is plagued with unphysical behavior, and, as we have 
shown, it is preferable to use a simplified version where 
occupation numbers are frozen to their initial value. 

In summary, both TDHFB and TDHF+BCS with con- 
stant occupation numbers can eventually be used to de- 
scribe a physical system while TDHF+BCS with varying 
occupation numbers should be avoided. While TDHFB is 
expected to have a reacher dynamics, due to its simplic- 
ity, the second transport theory remains quite attractive. 
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